What is the lowest possible reheating temperature? 
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We study models in which the universe exits reheating at temperatures in the MeV regime. By 
combining light element abundance measurements with cosmic microwave background and large 
scale structure data we find a fairly robust lower limit on the reheating temperature of Trh ^ 4 
MeV at 95% C.L. However, if the heavy particle whose decay reheats the universe has a direct 
decay mode to neutrinos, there are some small islands left in parameter space where a reheating 
temperature as low as 1 MeV is allowed. The derived lower bound on the reheating temperature 
also leads to very stringent bounds on models with n large extra dimensions. For n = 2 the bound 
on the compactification scale is M > 2000 TeV, and for n = 3 it is 100 TeV. These are currently 
the strongest available bounds on such models. 
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I. INTRODUCTION 

The standard big bang model has been tested thor- 
oughly up to temperatures around 1 MeV where big bang 
nucleosynthesis occurred. At much higher temperatures 
the universe is assumed to have undergone inflation, dur- 
ing which the primordial density perturbations are pro- 
duced. 

Towards the end of inflation the inflaton potential 
steepens so that slow roll is violated, and the universe 
enters the reheating phase. During this phase all parti- 
cles which are kinematically allowed are produced, either 
by direct decay or from the thermal bath produced by 
the inflaton decay. 

Finally the universe enters the radiation dominated 
phase at a temperature Trh, which is a function of the 
inflaton decay rate. The only certain bound on this re- 
heating temperature comes from big bang nucleosynthe- 
sis, and has in several previous studies been found to be 
around 1 MeV 

It should be noted that even if the reheating tempera- 
ture after inflation is much higher there can still be sub- 
sequent "reheating" phases, in the sense that reheating is 
defined to be a period where the energy density is domi- 
nated by an unstable non-relativistic particle species. In 
standard reheating this is the inflaton, but in supersym- 
mctric models it could for instance be the gravitino. 

In the present paper we update previous calculations 
of this reheating phenomenon, using data from cosmic 
microwave background and large scale structure observa- 
tions. Furthermore we extend the analysis to include the 
possibility of having a direct decay mode of the heavy 
particle into light neutrinos. If the heavy particle is 



a scalar this decay is normally suppressed by a factor 
(mu/m^)'^ because of the necessary helicity flip. How- 
ever, the heavy particle could either be a non-scalar parti- 
cle, or it could be a pseudo-scalar like the majoron which 
couples only to neutrinos. Even though such models are 
slightly contrived it is of interest to study whether the 
temperature bound on reheating is signiflcantly affected 
by the possibility of direct decay into neutrinos. 

In section II we discuss the set of Boltzmann equations 
necessary to follow the evolution of all particle species. 
In section III present results of the numerical solution 
of these equation, and in section IV we compare model 
predictions with observational data. Finally, section V 
is a review of other astrophysical constraints on heavy, 
decaying particles, and section VI contains a discussion. 



II. BOLTZMANN EQUATIONS 

We follow the evolution of all particles by solving the 
Boltzmann equation for each species 



df „ df 



(1) 



where Ccou is the collision operator describing elastic and 
inelastic collisions. 



A. Neutrinos 

Neutrinos interact with the electromagnetic plasma via 
weak interactions. A comprehensive treatment of this can 
for instance be found in Ref. 0|. The collision integrals 
can be written as 
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where 5*1-^112^34, i is the spin-summed and averaged ma- 
trix element including the symmetry factor S = 1/2 if 
there are identical particles in initial or final states. The 
phase-space factor is A(/i, /a, /3, fi) = fshC^ - /i)(l - 
/2)-/i/2(l-/3)(l-/4). 

This collision integral can be reduced to 2 dimensions 
using the method developed in Ref. However, if Pauli 
blocking and interactions involving only neutrinos are ne- 
glected the integrals can in fact be reduced to 1 dimen- 
sion, as described in Ref. 0. In the following we use 
this method. The quantitative error resulting from this 
is quite small. 

In addition to standard weak interactions we allow for 
a direct decay of 4> to neutrinos, (/> — > vv. If is non- 
relativistic then each neutrino is born with momentum 
and in this case the collision integral is 



:T^n^S{pi, ~ (3) 



where b^- is the branching ratio into neutrino species i, 
and is the decay rate of the heavy particle. For sim- 
plicity we assume equal branching ratios into all neutrino 
species. Even is this is not the case the neutrino distribu- 
tion functions will be almost equilibrated by oscillations 
0. This means that ~ — b,^^ ~b^/3. 

Note that if one assumes that neutrinos are in kinetic 
equilibrium so that they can be described by a single tem- 
perature it is in fact possible to solve the Boltzmann 
equation semi-analytically 5j. However, this is a very 
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poor approximation for the case when there is a direct 
decay mode (/) — > vv. 



We assume the heavy particle to be completely non- 
relativistic. If that is the case then the Boltzmann equa- 
tion can be integrated to give the following equation for 
the evolution of the energy density 



(4) 



i.e. there are no inverse decays. This is a good approxi- 
mation for all the cases covered in the present work. 

We only work with masses which are low enough that 
there are no hadronic decay channels open. This of course 
severely restricts the possible models. However, if there 
is a hadronic branching ratio then the minimum allowed 
reheating temperature increases dramatically , and we 
are investigating what the lowest possible reheating tem- 
perature is. 

C. Electromagnetic plasma 

The evolution of the photon temperature can then be 
found from the equation of energy conservation 



dt 



= -3H{pt + Pt), 



(5) 



where pT and Pt are the total energy density and the to- 
tal pressure respectively. This equation can be rewritten 
as an evolution equation for Tj 



dT-y 

~dr 



-(1 - b^)p^r^ + 4Hp-, + 3H{pe + Pe) + 4Hp^ + dp^/dt 



dp^/dT^ 

r 



dpe/BT^ 



(6) 



D. Scale factor 



E. Initial conditions 



Finally we solve the Friedmann equation to find the 
scale factor as a function of time 



H 



a / SnGpT 
a 



(7) 



Altogether we solve Eq. ^ together with Eq. |^ for 
each neutrino species, Eq. |0J for 0, and Eq. © for the 
photon temperature, to obtain a(t), T^{t), Pcj,{t), and 



Following convention we define the reheating temper- 
ature of the universe to be when 



= 3i7(rRH) 



(8) 



To a reasonable approximation the universe is radiation 
dominated at this point so that 



H = 



g*7r 
90 



2\ 1/2 ^2 

T 



RH 



(9) 



where Mp\ = 2.4 x 10^* GcV is the reduced Planck mass 
and is the number of degrees of freedom. 

This means that there is a one to one correspondence 
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between and Trh, 

TR,H,McV-0.7r'/?, (10) 

where = 10.75 has been used. Note that the con- 
stant of proportionahty is somewhat arbitrary (although 
it should always be of order 1), and just gives a rough 
idea about the thermal temperature when the universe 
enters the standard radiation dominated phase. 

As long as the initial time is set so that ti <^ ^(Trh) 
and Tmax ^ Tu,^, where Tmax is the maximum temper- 
ature reached by the plasma after time U and Td.i/ is 
the neutrino decoupling temperature then the final out- 
come is independent of initial conditions. The universe 
starts out being strongly matter dominated and the fi- 
nal neutrino energy density, as well as the light element 
abundances depend only on T^, m^, and b^. The initial 
time is found from the Friedmann equation by assuming 
complete domination of <j) so that U — |[87rGp0_i/3]~^/^. 

F. Nucleosynthesis 

One of the main observables from the epoch around 
neutrino decoupling is the abundance of light elements, 
mainly helium and deuterium. In order to calculate these 
abundances we have modified the Kawano nucleosynthe- 
sis code . First is has been modified to incorporate the 
modified temperature evolution, and second the subrou- 
tines used to calculate weak interaction rates for n <-> p 
have been modified to incorporate the full numerical elec- 
tron neutrino distribution coming from the solution of the 
coupled Boltzmann equations. 

This allows us to calculate the abundance of ^He and 
D for the various models. 



III. NUMERICAL RESULTS 
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FIG. 1: The effective number of neutrino species as a function 
of when there is no direct decay into neutrinos, bv = 0. 

for = 6.4 s^^ for two different initial times, U = 
1.8 X 10^^ s and ti ~ 8.8 x 10~^ s. In both cases we 
assume an initial photon temperature of 2.3 MeV (we 
could equally well have chosen an initial temperature of 
0). While the maximum temperature reached is clearly 
dependent on ti , T-y and quickly become indistinguish- 
able, and as long as the temperature where this happens 
is greater than the neutrino decoupling temperature all 
final results are independent of ti. Furthermore, as ex- 
pected [3|, the photon temperature scales as oc t~^^^ 
during the matter dominated period and shifts to the 
usual oc t~^^'^ once the universe becomes radiation 
dominated (except for a small deviation due to heating 
by e"^e~ annihilation). 



We have solved the set of coupled Boltzmann equations 
for all species for the free parameters, m<^, F^, and 6^. 

The main output from this is the relativistic energy 
density in neutrinos, parameterized in units of the energy 
density of a standard model neutrino, p^^ , 



A. b^=0 

lib,y — then the equations become independent of 
and this case has already been covered in Ref. 0. We 
present this as our first case in order to compare results 
with those of . Fig. Q] shows the effective number of 
neutrino species, N^, after complete decay of </>. This 
figure is identical to Fig. 4 in Ref. 0. 

We also test whether our results are independent of 
initial conditions. In Fig. |21 we show Ty{t) and p,p{t) 



B. 6, / 

Next we cover the case when b^, ^ 0. This is much more 
complicated to solve numerically because of the presence 
of the delta function d{p^ — and the fact that the 

solution now depends on both b^ and m^. In Fig. |31we 
show Ni, for different values of F^ and 5^ . 

From this figure is clear that when bi, is small the ef- 
fective number of neutrino species becomes independent 
of and increasing with F^, with N^, — > 3 for F^ ^ oo. 

For the opposite case when bi, = I (only decay to neu- 
trinos) the situation is the opposite. When F — > cxd the 
limiting value is again N^, = 3. This corresponds to the 
case when (j) decays into neutrinos, but the effective neu- 
trino temperature after complete (f) decay is higher than 
Td. 

When F ^ the effective number of neutrino species 
goes to infinity. This corresponds to the case when de- 
cays so slowly that the produced neutrinos never equili- 
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FIG. 2: and as functions of time for = 6.4 s~^, 
— and two different initial times. The full line is for 
U = 8.8 X 10~^ s, whereas the dashed is for U = 1.8 x 10^"* s. 



distribution functions. In Fig. ^ which shows 6^ = 1, 
— 6.4 s^^, and — 120 MeV, it can be seen that 
the distribution function is higher than thermal at high 
energies because of (j) decay. However, there are fewer low 
energy neutrinos because of the inefficient production via 
e+e^ annihilation. 

Conversely, in Fig.|Sl which shows b^, = 1, = 50 s~^, 
and = 120 MeV, it can be seen that the decay rate 
is high enough that neutrinos equilibrate with the elec- 
tromagnetic plasma, except for a small deviation around 
Pi, — m0/2. This subsequently leads to ~ 3 after 
complete (j> decay. 



IV. COMPARISON WITH DATA 

In order to constrain the parameters b^, F^, and 
we compare the predicted values of N,y, '*IIe, and D with 
the observationally determined values. In addition to 
the parameters directly related to (j) the nucleosynthesis 
outcome depends crucially on the baryon density, rj = 

Taken at face value the recent CMB data from the 
WMAP satellite constrain rj tightly. However, it has been 
shown that there is a significant correlation between 77 
and Ni, in the CMB data. This means that it is not 
possible to take CMB constraint on rj directly and apply 
it to the nucleosynthesis calculations. Rather a full CMB 
likelihood analysis for and r] must be carried out. This 
can then be combined with the nucleosynthesis likelihood 
analysis for 6^, F^, m^, and rj. 

First the following subsection covers the current obser- 
vational status, then the next covers the constraints on 
decay parameters which can be obtained. 



brate with the electromagnetic plasma, leaving only neu- 
trinos. 

However, there is a large intermediate region where 
Ni, < 3, even for by = I. The reason for this unexpected 
feature can be explained as follows: When high energy 
neutrinos {E ^ T) are produced by direct (f) decay they 
have a very high annihilation cross section to e"'"e~, be- 
cause the cross section goes as E"^ . However, the pro- 
duced electrons and positrons are immediately converted 
into a sea of low energy e"*", e~, and 7 because of electro- 
magnetic interactions. This means that the production 
rate of neutrinos is much lower. In the case where the 
reheating temperature is very high this does not matter 
because the universe still has time to thermalize com- 
pletely after (j) decay. However, if Trh ^ 7b, this is not 
possible and the result is that < 3 because of the 
very efficient conversion of neutrinos into e~^e~ . Notice 
also that this effect becomes less pronounced when 
decreases because neutrinos are born with energies closer 
to 3T, and the mismatch between forward and backward 
rates becomes smaller. 

In Figs. ^ and |S1 this effect can be seen directly on the 



A. Observational data 

1. Light element abundances 

The primordial helium abundance has been derived by 
two independent groups. Fields and Olive find the 
value 

= 0.238 ±0.002 ± 0.005, (12) 

whereas Izotov and Thuan find 

Fp = 0.244 ± 0.002 ± 0.005 (13) 

Because of this inconsistency we blow up the error bars 
on Yp and use the value 

yp = 0.238 ±0.015, (14) 

which encompasses the allowed regions of both observa- 
tional determinations. 

The most recent determination of the primordial deu- 
terium abundance has yielded the value [3 

D/H= (2.78 ±0.29) X 10"^ (15) 
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FIG. 3: Contour plot of Nu for different and T^. The top left plot is for bv = 0.1, the top right for bv = 0.5, the bottom 
left for bv — 0.9, and the bottom right for 6^ — 1.0. 



A determination of the primordial lithium abundance 
has also been performed by several groups. However, this 
measurement is prone to large systematics and we refrain 
from using it here. 



2. Cosmic microwave background 

The CMB temperature fluctuations are conveniently 
described in terms of the spherical harmonics power spec- 
trum 

Ci EE (|az„p), (16) 

where 

AT" 

— (0,0) =^a,™yi™(0,0). (17) 

Im 

Since Thomson scattering polarizes light there are ad- 
ditional power spectra coming from the polarization 



anisotropics. The polarization can be divided into a curl- 
free (E) and a curl (B) component, yielding four indepen- 
dent power spectra: Ct.i, Ce.i, Cb,i and the temperature 
-E-polarization cross-correlation Cte,i- 

The WMAP experiment have reported data on Ct,i 
and Cte,i, as described in Ref. pHlllLll^ 

We have performed the likelihood analysis using the 
prescription given by the WMAP collaboration which in- 
cludes the correlation between different Ci's [TollTllIT^ . 
Foreground contamination has already been subtracted 
from their published data. 



3. Large scale structure 

The 2dF Galaxy Redshift Survey (2dFGRS) [3 has 
measured the redshifts of more than 230 000 galaxies 
with a median redshift of « 0.11. An initial estimate 
of the convolved, redshift-space power spectrum of the 
2dFGRS has been determined [lj| for a sample of 160 
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FIG. 4: The distribution function for for different values 
of T-y when V = 6.4 s"\ K = 1, and = 120 MeV. The 
dotted line is for = 2.18 MeV, the dashed for = 0.42 
MeV, the long-dashed for = 0.19 MeV, and the full line 
for T-y = 0.01 MeV. The full grey (red) line is an equilibrium 
distribution with T,j = T-y. 



spectrum may be biased with respect to the matter power 
spectrum, i.e. light does not trace mass exactly at all 
scales. This is often parametrised by introducing a bias 
factor 



(18) 



where Pg{k) is the power spectrum of the galaxies, and 
Pmik) is the matter power spectrum. However, we re- 
strict our analysis of the 2dFGRS power spectrum to 
scales k < 0.15 liMpc"^ where the power spectrum 
is well described by linear theory. On these scales, 
two different analyses have demonstrated that the 2dF- 
GRS power spectrum is consistent with linear, scale- 
independent bias Thus, the shape of the galaxy 
power spectrum can be used straightforwardly to con- 
strain the shape of the matter power spectrum. 

The only parameters which affect CMB and structure 
formation are the baryon density, rj, and the relativistic 
energy density at late times, parameterized by Ni, [TtIIi^ 
(see also [19," 20, 21, 22, 23, 24, 25, 2^). It is therefore 
relatively straightforward to perform the CMB+LSS like- 
lihood analysis. 
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FIG. 5: The distribution function for for different values 
of when T = 50 s"\ K = 1, and = 120 MeV. The 
dotted line is for T-y = 7.7 MeV, the dashed for T-y = 0.93 
MeV, the long-dashed for T-y = 0.23 MeV, and the full line 
for T-y = 0.01 MeV. The full grey (red) line is an equilibrium 
distribution with T^ = T-y. 



000 redshifts. On scales 0.02 < fc < 0.15ft. Mpc"^ the 
data are robust and the shape of the power spectrum is 
not affected by redshift-space or nonlinear effects, though 
the amplitude is increased by redshift-space distortions. 
A potential complication is the fact that the galaxy power 



B. Likelihood analysis 

Nucleosynthesis is affected both by the expansion rate 
around T ^ 0.1 — 1 MeV, and by the electron neutrino dis- 
tribution function. The reason is that electron neutrino 
enter directly in the weak reactions which interconvert 
protons and neutrons. 

The specific neutrino distributions are therefore found 
as functions of temperature and used in a modified ver- 
sion of the Kawano BBN code @. This is then used 
to calculate primordial abundances of deuterium and he- 
lium. 

For calculating the theoretical CMB and matter power 
spectra we use the publicly available CMBFAST package 
j27l |. As the set of cosmological parameters we choose 
rim, the matter density, O^, the baryon density, Hq, the 
Hubble parameter, r, the optical depth to reionization, 
Q, the normalization of the CMB power spectrum, b, 
the bias parameter, and the effective number of neutrino 
species N^, found from the solution of the Boltzmann 
equations. We assume neutrinos to be almost massless. 
We restrict the analysis to geometrically flat models ilm+ 

For each individual model we calculate ij^ the fol- 
lowing way: Given a theoretical CMB spectrum the 
of the WMAP data is calculated using the method de- 
scribed in Ref. 0. With regards to the 2dF data we 
use the data points and window functions from Ref. p3 | 
(http://www.hcp.upenn.edu/^max/2df.html). 68% and 
95% confidence levels from the data are calculated from 
Ax^ = 2.31 and 6.17 respectively. 
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freeze-out there are many more protons than neutrons. 
When ^ m„ — nip the weak absorption cross section 
is equal on protons and neutrons. This means that addi- 
tional neutrinos at high energies will have the net effect of 
converting protons into neutrons, so that in the end more 
helium is produced. Note that this is the opposite effect 
of just increasing the weak interaction rates, in which 
case less helium would be produced. The phenomenon is 
quite similar to what happens if (f> has a hadronic decay 
channel. In that case pions and kaons will be produced, 
which subsequently convert protons to neutrons and lead 
to overproduction of helium. 

In Fig. we show 68%, 95%, and 99.99% confidence 
exclusion plots for and m^, marginalized over rj. 

Both when b^, is small and when = I the bound on 
7rh becomes independent of m^. In both cases the 95% 
bound is Trh ^ 4 MeV. 

However, there is an intermediate regime for bi, which 
allows for much lower values of Trh. The reason for 
this can be seen directly from Fig. O i.e. there is an 
intermediate range where can be kept close to 3, even 
for low F^. However, for large masses (which is of course 
by far the most likely) there is no allowed region. The 
reason is the one given in the previous section: More high 
energy neutrinos will produce more helium, and this in 
turn will conflict with observations. 

The final outcome is that for almost all values of 
and bi, there is a robust lower bound on Trh which is 
around 4 MeV. However there is a small region where 
b^ ~ 0.9, < 40 MeV where a reheating temperature 
as low as roughly 1 MeV is allowed. 



V. OTHER CONSTRAINTS 



FIG. 6: 68% and 95% confidence exclusion plot of the param- 
eters rjio = 10^" X ri and F^ for the case when b^, = 0. 



1. b„ = 

In Fig. we show 68% and 95% exclusion limits for 
T] and F^ from BBN, CMB, and LSS. The top panel 
for BBN only is very similar to Fig. 8 in KKS, except 
that we use slightly different bounds on light element 
abundances. From BBN alone the 95% bound on Trh is 
roughly 0.6 MeV. However this bound is achieved for rel- 
atively low 77, whereas CMB-I-LSS strongly prefer a high 
value of 77. Therefore combining the BBN and CMB-I-LSS 
constraints removes the low Trh region and increases the 
lower bound to 3.9 MeV. 



2. 6, / 

Apart from the fact that N^, depends on 6,^ there is a 
second effect which is just a important. When b^, 
there are more high energy neutrinos. Around weak 



If (/> is a scalar then the decay rate to neutrinos is 
normally suppressed by a factor because of the nec- 
essary helicity flip. Therefore the simplest assumption 
is that (f) has no branching into neutrinos. If for in- 
stance the heavy particle is a pseudo-scalar like the axion, 
then there is an upper bound on the coupling to photons 
(113, EH), (707 < 0.6x10-1" GcV-i forTO0 < 30keV. For 
higher masses the bound is significantly weaker. How- 
ever, even if this bound is used together with the decay 
width F0^27 = 30-Y™0/647r then we find that 

r0^27 ^ SOmj^ loMoV (19) 

which is easily satisfies for the parameter space we are 
considering. 

On the other hand, if </> is a particle like the majoron 
which couples only to neutrinos then the decay width is 

m 

2 

F^^.p = ^^^^ ~ 3 X IQi^ m^.McV s'^ (20) 

The bound on the dimensionless coupling constant comes 
from BBN as well as supernova considerations and is of 
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FIG. 7: 68%, 95%, and 99.99% confidence exclusion plot of the parameters and using all available data 

(CMB+LSS+BBN). The top left plot is for K = 0.1, the top right for = 0.5, the bottom left for K = 0.9, and the 
bottom right for 6^ = 1.0. 



order 10^^ — 10^'"' for majorons in the MeV mass range 
|8lLl8.'jl |. For more massive majorons the bound weakens. 
Again it is clear that the decay parameters which we 
consider here are not excluded by any other astrophysical 
or experimental data. 

The final conclusion is that heavy, decaying particles 
such as the ones considered here cannot be directly ex- 
cluded by any current data. Furthermore a branching 
ratio into neutrinos can be anywhere from to 1. 



VI. CONCLUSION 

We have carefully calculated constraints on models 
with extremely low reheating temperature, where a mas- 
sive particle decays around T ^ 1 MeV. By combining 
constraints on light element abundances with constraints 
on r] and N^, from CMB and large scale structure we 



derived a fairly robust limit of 

Trh > 4 MeV. (21) 

This bound is a significant improvement over the pre- 
vious bound of Trh ^ 0.7 MeV, calculated from BEN 
alone. It is interesting that the lower bound is signifi- 
cantly higher than the n <-> p conversion freeze-out tem- 
perature, T ~ 0.8 MeV, end even higher than the neu- 
trino decoupling temperature Td ~ 2 MeV. This shows 
that even small residual effects can be measured with 
present observational data. 

Models with reheating temperature in the MeV regime 
are in general difficult to reconcile with such features 
as baryogenesis. However, in models with large extra 
dimensions a low reheating temperature is essential in 
order to avoid overproduction of massive Kaluza-Klein 
graviton states. This means that we can use our present 
bound to derive limits on the compactification scale in 
such models. For the case of two extra dimensions the 
bound is M > 2000 TeV and for n = 3 it is M > 100 
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TeV. This bound is somewhat stronger than the bound 
coming from considerations of neutron star coohng and 
gamma ray emission. 



darriaga [27| , as well as the nucleosynthesis code written 
by Lawrence Kawano 0. I wish to thank P. Serpico for 
comments. 
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